The distribution of the ratio of consecutive level spacings in random matrix ensembles 
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We derive expressions for the probability distribution of the ratio of two consecutive level spacings for the 
classical ensembles of random matrices. This ratio distribution was recently introduced to study spectral prop- 
erties of many-body problems, as, contrary to the standard level spacing distributions, it does not depend on the 
local density of states. Our Wigner-like surmises are shown to be very accurate when compared to numerics 
and exact calculations in the large matrix size limit. Quantitative improvements are found through a polynomial 
expansion. Examples from a quantum many-body lattice model and from zeros of the Riemann zeta function 
are presented. 
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Random matrix theory (RMT) was introduced half a cen- 
tury ago in order to describe statistical properties of energy 
levels of complex atomic nuclei HI . Since then, it has proven 
to be very useful in a great variety of different fields Q2J, |3| . 

In quantum chaos fl4[], RMT accurately accounts for the 
spectral statistics of systems whose classical counterpart is 
chaotic. While for quantum Hamiltonians which classical 
counterpart is integrable, the Berry-Tabor conjecture [5] states 
that their level statistics follows a Poisson law, Bohigas, Gi- 
annoni and Schmit conjectured [6] that the case of quantum 
Hamiltonians with chaotic classical dynamics must fall into 
one of the three classical ensembles of RMT. These three 
ensembles correspond to Hermitian random matrices whose 
entries are independently distributed respectively real (GOE), 
complex (GUE) or quaternionic (GSE) random variables (see 
H for details). 

Universality of RMT means that random matrix ensembles 
describe energy levels of real systems at a statistical level, and 
only in a local energy window when the mean level density 
is set to unity. Different models may and do have very differ- 
ent level densities and to compare usual spectral correlation 
functions like the nearest-neighbor spacing distribution one 
has to perform a transformation called unfolding [Q] Q • The 
unfolding procedure consists in changing variables from the 
true levels, e„, to new ones, e„ = AT(e„), where A/"(e) is the 
mean number of levels less than e, obtained either by smooth- 
ing over many realizations in the case of disordered systems, 
or by local smoothing over an energy window large compared 
to the level spacing but small compared to variations of J\F(e). 
The unfolded spectrum has automatically a mean level spac- 
ing equal to one, and its statistical properties can thus be di- 
rectly compared with those of RMT. When a functional form 
of N is known (as for billiards), or when large enough statis- 
tics is available, the unfolding is straightforward and easily 
implemented. 

The situation is different for many-body problems, where 
AT(e) increases as a stretched exponential function of energy 
10] with, in general, unknown lower-order terms, and where it 
is difficult to calculate a large number of realizations because 
of an exponential increase of the Hilbert space dimension with 
the number of particles. In order to circumvent these diffi- 



culties which greatly diminish the precision of statistical tests 
in systems with a large number of particles, Oganesyan and 
Huse |Ql proposed a new quantity defined as follows. Let e n 
be an ordered set of energy levels and s n — e n+ i — e n the 
nearest-neighbor spacings. Oganesyan and Huse considered 
the distribution of the ratios f n defined by 



min(s n ,s n _i) . / 1 

r n = ; r = mm r„, — 

max s n ,s„_i V r„ 
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This quantity has the advantage that it requires no unfolding 
since ratios of consecutive level spacings are independent of 
the local density of states. Such a distribution thus allows a 
more transparent comparison with experiments than the tradi- 
tional level spacing distribution. For this reason, many recent 
works use this quantity in different contexts of many-body 
systems. As an example let us mention quantum quenches, 
where the tools of RMT and quantum chaos were used as a 
phenomenological approach to quantify the distance from in- 



, and also to investigate 
T3l . In these papers 



tegrability on finite size lattices [9-1 1 
numerically many-body localization 
the distribution of consecutive level spacing ratios P(f) was 
shown to yield more precise results than the usual spacing dis- 
tribution P(s). 

Although the distribution P(r) plays a more and more im- 
portant role in the interpretation of numerical data in quantum 
many -body Hamiltonians, only numerical estimates of it exist, 
and they are restricted to the GOE ensemble. RMT predictions 
for P{r) are lacking. Such predictions are essential, both for 
understanding its shape for the three RMT ensembles, and for 
providing accurate estimates with simple formulas that could 
be used as an efficient tool. 

This letter fills this gap by providing several important re- 
sults on P(r). First, we compute Wigner-like surmises for all 
three classical RMT ensembles, which already provide sim- 
ple analytical formulae in very good agreement with exact 
numerics and analytical expressions in the large matrix size 
limit. Second, the remaining small differences are shown to 
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be well fitted to numerical precision by a rather simple poly- 
nomial expansion. Results are then applied to examples on a 
quantum many -body Hamiltonian and to zeros of the Riemann 
zeta function. 

The ratio of consecutive level spacings distribution - In- 
stead of the quantity (Q]), we find it more natural to consider 
directly the ratio of two consecutive level spacings (fJJ and 
its probability distribution P(r). Indeed, let p{e\, e^,e^) be 
the probability density of three consecutive levels with ei < 
e 2 < e,3. Assuming translation invariance, p(ei, e2, 63) — 
P(si, 82) where s l = e i+ x - e». Then 

P(r) = J P(s 1 ,s 2 )S (r - ^\ d Sl ds 2 

pOO 

= / P(rS 2 ,S 2 )s 2 ds 2 . 

Jo 
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It is physically natural and can be proved analytically that for 
all classical RMT ensembles in the bulk of the spectrum (as 
well as for Poisson variables) the function P(si, s 2 ) is sym- 
metric, that is, P(si, s 2 ) — P(s 2 > si). This left-right symme- 
try implies then that the distributions of r n and 1 / r„ are the 
same, so that P(r) satisfies the following functional equation 



p (r) = ±P (l 



(4) 



Whenever (O holds, it is equivalent to consider the whole dis- 
tribution P(r) or to restrict the study to the support [0, 1] by 
considering the variable f defined in (Q3, as was done in fl. 
Here we concentrate on the whole distribution P(r); since 
P(f) — 2P(r)0(l — r), our results can easily be translated to 
the restricted distribution. The integrable (Poisson) case triv- 
ially yields P(r) = 1/(1 + r) 2 . We now address the behavior 
of P(r) for RMT ensembles. 

Wigner-like surmise - For Gaussian ensembles, the joint 
probability distribution of N eigenvalues is given by yfl 



p(ei, . . .,e N ) 
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where Cp n is a known normalization constant and f3 is the 
Dyson index equal to 1 (GOE), 2 (GUE) or 4 (GSE). The ex- 
act calculation of P(r) via Eq. (0 requires the calculation 
of P(si,s 2 ). Though this calculation is possible from <j3j 
(as shown at the end of this Letter), it ultimately requires the 
use of numerical methods and is not transparent. Exactly the 
same problem appears in the calculation of the usual nearest- 
neighbor spacing distribution, P(s), which is the probability 
that the distance between two consecutive levels is s. Rather 
than cumbersome exact calculations, Wigner derived a simple 
approximate expression for P(s), 



P w (s) = ap 



(6) 



with some explicitly known normalization constants ap and 
bp J2l- This formula, called the Wigner surmise, corresponds 



to the exact result for 2x2 matrices, and is in very good 
agreement with the exact large- N expressions 111 311 . 

In a similar spirit, we obtain a formula for the ratio distri- 
bution of two consecutive spacings by performing the exact 
calculation for 3x3 matrices, starting from the joint dis- 
tribution (|5j for three eigenvalues ex, e 2 , e^. If for instance 
ei < e 2 < e3, the ratio r is given by (e 3 — e 2 )/(e2 — ei). 
Consequently, the distribution P(r) in the 3x3 case is pro- 
portional to 
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After the change of variables x = e 2 — ex, y = e 3 — e 2 , the 
integration over e 2 is trivial and the remaining integrals read 



dxdy6(rx - y)x f3+1 y fi {x + y f e -^ x2+y2)+ ^ x ~ v)2 . 



After performing the integrals, the surmise takes the simple 
form 



Pw(r) 
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Zp (l + r + r 2 ) 1+ i^ 



(7) 



with Zp the normalization constant (see values in Tabled. 

One can check that this result satisfies the symmetry (0). 
The distribution Pyv (r) has the same level repulsion at small r 
than P(s), namely Pvv (r) ~ r' 3 , while for large r the asymp- 
totic behavior is Pvv{ r ) ~ r~( 2+/3 \ contrary to the fast ex- 
ponential decay of P(s). This surmise also yields an analytic 
expression for the mean-values (r)w an d (f)w widely used 
in the literature as a measure of chaoticity (see Table |I]for the 
exact values). 

Comparison with numerics and polynomial fit - We now 
investigate the accuracy of the surmise (|7) with respect to nu- 
merical calculations for large matrix sizes. As illustrated in 
Fig. Q] the surmise is almost indistinguishable from numer- 
ics and can thus be used for practical purposes as a reference 
to discriminate between regular and chaotic dynamics. The 
absolute difference 6P(r) — P num (r) — Pwir) between nu- 
merics and the surmise d7} is plotted in Fig. [2] for the three 
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TABLE I. Values of useful constants and averages (r) and (f }. Aver- 
ages (.) w are calculated from Eq. ((7}, and (.)g t from data in Fig.[T] 
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FIG. 1. (Color online) Distribution of the ratio of consecutive level 
spacings P(r) for Poisson and RMT ensembles: full lines are the 
surmise Eq. ((7}, points are numerical results obtained by diagonal- 
izing matrices of size N = 1000 with Gaussian distributed entries, 
averaged over 10 histograms. Inset: the distribution P(r). 



ensembles, and has a maximum relative deviation of about 
5%, similar to the Wigner surmise for P(s) lfl3ll . 

In order to go beyond the surmise (|7}, we propose a sim- 
ple expression which perfectly fits this remaining difference 
5P(r) within our computational accuracy. In order to fulfill 
Eq. (|4j, and assuming that P(r) for large N and Pw(r) have 
the same asymptotic behavior for small and large r, a reason- 
able ansatz is the following expansion 



M%t(r) 



C 



1 



-0 
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-09+1) 



(8) 

where cp is easily calculated from the normalization condition 
L 5P(r)dr = (see Table J] for the exact value). Thus the 
large- N expression for P(r) can be fitted by the expression 
P(r) = Pw(r) + 5Pflt(r) with only one fitting parameter, 
which is the overall magnitude C of the discrepancy. The 
best fit C can be found in Table J] The corresponding curves 
are shown in Fig. [2] Thanks to these very good fits, one can 
quickly infer accurate predictions for (r) and (f) and any 
average weighted by P{r) (see Tabled). 

Large-N calculation - We now turn to the exact calculation 
of P(r) for GUE (i.e. /3 = 2) in the limit N -> oo, follow- 
ing a path similar to the derivation of the exact level spacing 
distribution P(s). 

Our starting point is Eq. 5.4.29 of Ref. [2. From that equa- 
tion, one can check that the probability p(— t, y, t) of having 
three consecutive levels at points — t, y, t can be rewritten as 

p(-t, y, t) = det(l - K) det[R(x, «) x ,*=-t,»,t], (9) 

where R(x, y) is the resolvent kernel, i. e. the kernel of the 
operator (1 — K)^ 1 K, and det(l — K) is the Fredholm deter- 
minant of K. Operator K is an integral operator whose action 



FIG. 2. (Color online) Difference 5P{r) = P mm (r) - P w (r) be- 
tween the numerics and the surmise ((T}. The fit function is given by 
Eq. (8). Green diamonds are results of exact calculations obtained 
from l|9} for GUE. 



is defined as 



(Kf)(x) 



with the kernel 



K(x,y) = 



K (x,y)f{y)dy 



sin7r(a; — y) 



tt(x - y) 



(10) 
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It is known (see e. g. Il4ll ) that for a kernel of this form the 
resolvent kernel can be written as 



R(x,y) = 



Q(x)P(y) - Q(y)P(x) 



x-y 



(12) 



with functions Q(x) and P(x) obeying integral equations 



Q(x) 



K(x,y)Q(y)dy = 



P{x)- J K(x,y)P(y)dy 



cos nx. 



(13) 



Function Q(x) and P{x) have many useful properties which 
allow to relate the calculation of spectral statistics for stan- 
dard RMT ensembles to solutions of Painleve equations (see 
e. g- IU4I1 and references therein). Though this approach is ele- 
gant, it still requires numerical resolution of Painleve V equa- 
tion for det(l — K) with subsequent solutions of linear equa- 
tions for Q(x) and P(x) whose coefficients are determined by 
that solution. 



We find it simpler to use the direct method proposed in J15 1 
for computing dct(l — K). It is based on a quadrature method 
for numerical evaluation of the integrals 



f{x)dx = y^^Wkf(xk) 
fc=i 



(14) 



appearing in the definition (TTOb of the integral operator K. 
Such a discretization allows to approximate the determinant 
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FIG. 3. (Color online) (a) P(r) — Poo (r) for GUE and various matrix 
sizes. Inset: constant C'n from the fit l[8} as a function of matrix size 
TV (solid line is a fit 1/N). (b) Density distributions for the overlap- 
ping ratio r£ 2) = (e n+2 — e„)/(e n+ i — e„_i) for Poisson variables 
and for the three classical RMT ensembles (same color code as in 

Fig.rj. 



of the integral operator as a finite m x m determinant 

det(l - K) « det (5. jk - K(x j ,x k )wk) 



(15) 



and functions Q(x) and P(x) defined in (TOl ) can be ob- 
tained by solving a linear system of m equations. As noted 
in fl5f l the method quickly converges. The result is presented 
in Fig. |2j where the Clenshaw-Curtis method with up to 60 
points of discretization has been used for the discretization 
(TBI . Figure|3](left) shows how the numerical results converge 
to the analytic large- N calculation. As mentioned previously, 
the fit P(r) = Pw(r) + 6Pm(r) works well for all N, with 
an overall TV-dependent constant CV in (©. This constant, 
which gives the amplitude of the departure from the Wigner- 
like surmise, asymptotically decreases as 1/N (see inset of 
Fig.©. 

Applications - To illustrate the above formalism, we inves- 
tigate the spectral properties of a quantum Ising chain of L 
spins— \ with periodic boundary conditions in transverse field 
A and longitudinal field a. The Hamiltonian is given by 



H = — 



L 

£ 

n=l 



'n u n+l 



'L+l 



(16) 

where ct^ 2 are the Pauli matrices at site n. This model re- 
cently attracted attention due to its experimental realization in 
cobalt niobate ferromagnet I1 1611 . The Hamiltonian ([TBI com- 
mutes with the operator T which translates the state by one 
lattice spacing and obeys T L = 1. Consequently, H takes a 
block diagonal form in the basis of eigenstates of T, and one 
has to consider separately each sector of symmetry. The result 
for one sector is illustrated in Fig. |4] Other symmetry sectors 
give similar results. As expected, P(r) agrees well with the 
GOE prediction © with (3 = 1. 

Another example of application is to look at non-trivial ze- 
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FIG. 4. (Color online) Histogram of the ratio of consecutive level 
spacings P(r). Black: Quantum Ising model in fields A = a = 
0.5 in sector of eigenvectors of T with eigenvalue uus (ojj = 
exp(2i7rj/L), j = 0, 1, . . . , L — 1), for L — 18 spins (dimen- 
sion of eigenspace = 14541). Violet: The same for zeros of Riemann 
zeta function up the critical line (10 4 levels starting from the 10 22 th 
zero, taken from fUjll ). Full lines correspond to the Wigner-like sur- 
mise Eq. Q with respectively 13 — 1 and f3 = 2. Inset: Difference 
between the numerics and these surmises. 



ros of the Riemann zeta function 



as) 



oo 1 

Y- 

n=l 



(17) 



It is well established that statistical properties of Riemann ze- 
ros are well described by the GUE distribution [17]. The prob - 
ability distribution of the ratio of two consecutive spacings of 
these zeros, presented in Fig. |4] is in a perfect agreement with 
GUE formula Q with (3 = 2. 

Conclusion - The investigation of spectral statistics in 
many-body problems with a large number of particles at- 
tracted wide attention in recent years. The absence of a well 
established expression for the mean density of states greatly 
diminishes the usefulness of standard correlation functions 
such as the nearest-neighbor spacing distribution. To avoid 
this problem, a new statistical tool has been proposed in il, 
namely the distribution of the ratio of two consecutive level 
spacings. 

The main result of the paper is the derivation of simple ap- 
proximative formulae for this distribution for classical RMT 
ensembles. The resulting Wigner-like surmises agree very 
well with direct numerical calculations. The difference be- 
tween the surmise and the exact calculations is small and can 
be fitted by a one-parameter polynomial formula with excel- 
lent accuracy. 

In the same spirit, several different ratios can be introduced 
which generalize the quantity ©. Analytic expressions and 
Wigner-like surmises can be derived in a similar way for the 
density distributions of these quantities, and will be discussed 
elsewhere. An example is given in Fig. (0 (right). All these 
distributions are universal in the sense that they apply with- 
out any unfolding or renormalization to spectra ranging from 
many-body systems to Riemann zeta function. 
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